!--------------------------------------------------------------------------------------------
! FileName:	gol_mpi.F90
! 
! Compile :	$ mpif90 -g -o gol_mpi_p2p.x -Wall -D_MPI_P2P gol_mpi.F90
!		$ mpif90 -g -o gol_mpi_bcast.x -Wall -D_MPI_BCAST gol_mpi.F90
!		$ mpif90 -g -o gol_mpi_data_derived.x -Wall -D_MPI_DERIVED_DATA gol_mpi.F90
!
! Run	  :	$ mpirun -np 4 -hostfile hosts ./gol_mpi_p2p.x
!
! Be careful:	1. GridX, GridY must be same size. (ex) 10x10, 20x20, 100x100
!		2. To do preprocessor, .F90 is required.
!		3. # of processor is 4, 9, 16, 25....f
!
! Tab size:	8
! Version :	1.12
!
! Updated by 	2011. 10. 14
!--------------------------------------------------------------------------------------------

PROGRAM gol_mpi_p2p
IMPLICIT NONE
INCLUDE 'mpif.h'
	INTEGER :: ROOT = 0
	INTEGER	nRank, nProcs, iErr, status(MPI_STATUS_SIZE)
	INTEGER, PARAMETER ::	nGridXSize=20, nGridYSize=20, nSteps=100, nSleep=3
	INTEGER nSqrtNprocs, nMyWidthX, nMyWidthY, nLeftOver, nQuotient, nRemainder
	INTEGER	nMyStartX, nMySTartY
	INTEGER	nTmpMyStartX, nTmpMyStartY, nTmpMyWidthX, nTmpMyWidthY
	INTEGER	i, j, k, l
	INTEGER	nLiveCell
	
	INTEGER, DIMENSION(:, :), ALLOCATABLE			:: gol_buf
	INTEGER, DIMENSION(:, :), ALLOCATABLE			:: tmp_gol_buf
	
#if defined (_MPI_DERIVED_DATA)
	INTEGER	column_new_type
#endif
	
	! MPI Start	
	CALL MPI_INIT(iErr)
	CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nProcs, iErr)
	CALL MPI_COMM_RANK(MPI_COMM_WORLD, nRank, iErr)
	
	! Let each processor decide what column(s) to run	
	IF (nProcs == 2) THEN
		nSqrtNprocs = 2
	ELSE
		nSqrtNprocs = SQRT(REAL(nProcs))
	ENDIF
	
	nMyWidthX = nGridXSize / nSqrtNprocs
	nMyWidthY = nMyWidthX
	nLeftOver = MOD(nGridXSize, nSqrtNprocs)
	
	nQuotient = nRank / nSqrtNprocs
	nRemainder = MOD(nRank, nSqrtNprocs)
	
	nMyStartX = nQuotient * nMyWidthY
	nMyStartY = nRemainder * nMyWidthY
	
	! Right edge
	IF (MOD(nRank+1, nSqrtNprocs) == 0) THEN
		nMyWidthY = nMyWidthY + nLeftOver
	ENDIF
	
	! Bottom edge
	IF (nProcs == 2) THEN
		nMyWidthX = nGridYSize
	ELSEIF ((nRank/nSqrtNprocs) == (nSqrtNprocs-1)) THEN
		nMyWidthX = nMyWidthX + nLeftOver
	ENDIF

	ALLOCATE(gol_buf(0:nGridXSize, 0:nGridYSize))
	
	! set to 0 
	DO i=0, nGridXSize-1
		DO j=0, nGridYSize-1
			gol_buf(i,j) = 0
		END DO
	END DO
	
	if (nRank == ROOT) THEN
		WRITE (*, *) 'nProcs = ', nProcs
		WRITE (*, *) 'nRank = ', nRank
		WRITE (*, *) 'nSqrtNprocs = ', nSqrtNprocs
		WRITE (*, *) 'nMyWidthX = ', nMyWidthX
		WRITE (*, *) 'nMyWidthY = ', nMyWidthY
		WRITE (*, *) 'nLeftOver = ', nLeftOver
	
		CALL getRandomValue(gol_buf, nGridXSize, nGridYSize)
		
		!DO i=0, nGridXSize-1
		!	DO j=0, nGridYSize-1
		!	PRINT *, i, j, gol_buf(i,j)
		!END DO
		
		CALL showGridState(gol_buf, nGridXSize, nGridYSize)
		WRITE(*, *)	
	ENDIF
	
	! Allocate memory buffer for holding own GOL values
	IF (nRank .NE. ROOT) THEN
		ALLOCATE(tmp_gol_buf(0:nMyWidthX, 0:nMyWidthY))
	END IF

	! -----------------------------------------------------------------------------------------------------
	! Main Loop : time steps
	! -----------------------------------------------------------------------------------------------------
	DO i=1, nSteps
		CALL MPI_BARRIER(MPI_COMM_WORLD, iErr)

		! Let all processor have same data		
#ifdef _MPI_P2P
		IF (nRank == ROOT)	WRITE (*, *) 'msg> #ifdef _MPI_P2P'
		
		IF (nRank == ROOT) THEN
			DO j=1, nProcs-1
				DO k=0, nGridXSize-1
					DO l=0, nGridYSize-1
						CALL MPI_SEND(gol_buf(k, l), 1, MPI_INTEGER, j, 10, MPI_COMM_WORLD, iErr)
					END DO
				END DO
			END DO
		ELSE
			DO j=0, nGridXSize-1
				DO k=0, nGridYSize-1
					CALL MPI_RECV(gol_buf(j, k), 1, MPI_INTEGER, ROOT, 10, MPI_COMM_WORLD, status, iErr)
				END DO
			END DO
		END IF
		
#elif defined (_MPI_BCAST)
		IF (nRank == ROOT)	WRITE (*, *) 'msg> #elif defined (_MPI_BCAST)'
		
		DO j=0, nGridYSize-1
			CALL MPI_BCAST(gol_buf(0, j), nGridXSize, MPI_INTEGER, ROOT, MPI_COMM_WORLD, iErr)
		END DO

#elif defined (_MPI_DERIVED_DATA)
		IF (nRank == ROOT)	WRITE (*, *) 'msg> #elif defined (_MPI_DERIVED_DATA)'
		
		CALL MPI_TYPE_CONTIGUOUS(nGridXSize, MPI_INTEGER, column_new_type, iErr)
		CALL MPI_TYPE_COMMIT(column_new_type, iErr)
		
		DO j=0, nGridYSize-1
			CALL MPI_BCAST(gol_buf(0, j), 1, column_new_type, ROOT, MPI_COMM_WORLD, iErr)
			!CALL MPI_BCAST(gol_buf(0, j), nGridXSize, MPI_INTEGER, ROOT, MPI_COMM_WORLD, iErr)
		END DO

		!CALL MPI_TYPE_FREE(column_new_type)
#endif
		! Do work at the my own grid cell
		CALL calEvolution(gol_buf, nGridXSize, nGridYSize, nMyStartX, nMyStartY, nMyWidthX, nMyWidthY)

		! Send my own do-work data to ROOT(0)
		IF (nRank .NE. ROOT) THEN
			DO j=nMyStartX, nMyStartX+nMyWidthX-1
				DO k=nMyStartY, nMyStartY+nMyWidthY-1
					tmp_gol_buf(j-nMyStartX, k-nMyStartY) = gol_buf(j, k)
				END DO
			END DO

			CALL MPI_SEND(nMyStartX, 1, MPI_INTEGER, ROOT, 20, MPI_COMM_WORLD, iErr)
			CALL MPI_SEND(nMyStartY, 1, MPI_INTEGER, ROOT, 21, MPI_COMM_WORLD, iErr)
			CALL MPI_SEND(nMyWidthX, 1, MPI_INTEGER, ROOT, 22, MPI_COMM_WORLD, iErr)
			CALL MPI_SEND(nMyWidthY, 1, MPI_INTEGER, ROOT, 23, MPI_COMM_WORLD, iErr)
			
#if defined (_MPI_P2P) || (_MPI_BCAST)
			IF (nRank == 1)		WRITE (*, *) 'msg> #if defined (_MPI_P2P) || (_MPI_BCAST)'	
			DO j=nMyStartX, nMyStartX+nMyWidthX-1
				DO k=nMyStartY, nMyStartY+nMyWidthY-1
					CALL MPI_SEND(tmp_gol_buf(j-nMyStartX, k-nMyStartY), 1, MPI_INTEGER, ROOT, 30, MPI_COMM_WORLD, iErr)
				END DO
			END DO
			
#elif defined (_MPI_DERIVED_DATA)
			IF (nRank == 1)		WRITE (*, *) 'msg> #elif defined (_MPI_DERIVED_DATA)'
			CALL MPI_TYPE_CONTIGUOUS(nMyWidthX, MPI_INTEGER, column_new_type, iErr)
			CALL MPI_TYPE_COMMIT(column_new_type, iErr)
			
			DO j=nMyStartY, nMyStartY+nMyWidthY-1
				CALL MPI_SEND(tmp_gol_buf(0, j-nMyStartY), 1, column_new_type, ROOT, 30, MPI_COMM_WORLD, iErr)
			END DO

			!CALL MPI_TYPE_FREE(column_new_type)
#endif
		ELSE IF (nRank == ROOT) THEN
			DO j=1, nProcs-1
				CALL MPI_RECV(nTmpMyStartX, 1, MPI_INTEGER, j, 20, MPI_COMM_WORLD, status, iErr)
				CALL MPI_RECV(nTmpMyStartY, 1, MPI_INTEGER, j, 21, MPI_COMM_WORLD, status, iErr)
				CALL MPI_RECV(nTmpMyWidthX, 1, MPI_INTEGER, j, 22, MPI_COMM_WORLD, status, iErr)
				CALL MPI_RECV(nTmpMyWidthY, 1, MPI_INTEGER, j, 23, MPI_COMM_WORLD, status, iErr)
				
				! WRITE(*, *) 'nRank=', nRank, 'nTmpMyStartX=', nTmpMyStartX, 'nTmpMyStartY=', nTmpMyStartY, 'nTmpMyWidthX=',&
				!	nTmpMyWidthX, 'nTmpMyWidthY=', nTmpMyWidthY
							
				ALLOCATE(tmp_gol_buf(0:nTmpMyWidthX, 0:nTmpMyWidthY))
				
#if defined (_MPI_P2P) || (_MPI_BCAST)
				IF (j == 1)	WRITE (*, *) 'msg> #if defined (_MPI_P2P) || (_MPI_BCAST)'					
				DO k=0, nTmpMyWidthX-1
					DO l=0, nTmpMyWidthY-1
						CALL MPI_RECV(tmp_gol_buf(k, l), 1, MPI_INTEGER, j, 30, MPI_COMM_WORLD, status, iErr)
					END DO
				END DO
#elif defined (_MPI_DERIVED_DATA)
				IF (j == 1)	WRITE (*, *) 'msg> #elif defined (_MPI_DERIVED_DATA)'
				CALL MPI_TYPE_CONTIGUOUS(nTmpMyWidthX, MPI_INTEGER, column_new_type, iErr)
				CALL MPI_TYPE_COMMIT(column_new_type, iErr)
				
				DO k=0, nTmpMyWidthY-1
					CALL MPI_RECV(tmp_gol_buf(0, k), 1, column_new_type, j, 30, MPI_COMM_WORLD, status, iErr)
				END DO
				
				!CALL MPI_TYPE_FREE(column_new_type)
#endif				

				DO k=nTmpMyStartX, nTmpMyStartX+nTmpMyWidthX-1
					DO l=nTmpMyStartY, nTmpMyStartY+nTmpMyWidthY-1
						gol_buf(k, l) = tmp_gol_buf(k-nTmpMyStartX, l-nTmpMyStartY)
					END DO
				END DO
			
				DEALLOCATE(tmp_gol_buf)
			END DO
		END IF

		CALL MPI_BARRIER(MPI_COMM_WORLD, iErr)
		
		IF (nRank == ROOT) THEN
			WRITE(*, *)
			WRITE(*, '(A, I6, A, I6)') 'msg> nStep = ', i, '/', nSteps
			WRITE(*, *)
			CALL showGridState(gol_buf, nGridXSize, nGridYSize)
	
			nLiveCell = 0
	
			DO j=0, nGridXSize-1
				DO k=0, nGridYSize-1
					nLiveCell = nLiveCell + gol_buf(j, k)
				END DO
			END DO
			
			WRITE(*, *)
			WRITE(*, '(A, I4)') 'msg> nLiveCell = ', nLiveCell
	
			WRITE(*, '(A, I2, A)') 'msg> Wait a ', nSleep, ' sec'
			CALL Sleep(nSleep)
		END IF
	END DO
	! -----------------------------------------------------------------------------------------------------
	! End of Main Loop : time steps
	! -----------------------------------------------------------------------------------------------------
	
	IF (nRank .NE. ROOT) THEN
		DEALLOCATE(tmp_gol_buf)
	END IF
	
	DEALLOCATE(gol_buf)

	! MPI End
	CALL MPI_FINALIZE(iErr)

!------------------------------------------------------------
CONTAINS
SUBROUTINE getRandomValue(gol_buf, nGridXSize, nGridySize)
INTEGER nGridXSize, nGridYSize
INTEGER, DIMENSION(0:nGridXSize, 0:nGridYSize)		:: gol_buf

	INTEGER	ii, jj
	
	!CALL RANDOM_SEED()
	
	CALL init_random_seed()
	
	DO ii=0, nGridXSize-1
		DO jj=0, nGridYSize-1
			gol_buf(ii, jj) = MOD(IRAND(), 2)	! Get random # (0 ~ 1) 
		END DO	
	END DO
END SUBROUTINE

SUBROUTINE showGridState(gol_buf, nGridXSize, nGridYSize)
INTEGER nGridXSize, nGridYSize
INTEGER, DIMENSION(0:nGridXSize, 0:nGridYSize)		:: gol_buf
	INTEGER	ii, jj
	
	WRITE(*, 400, advance='no')	'   '
	
	DO ii=0, nGridYSize-1
		WRITE(*, 100, advance='no')	ii
	END DO
	
	WRITE(*, *)				! new line
	WRITE(*, 400, advance='no')	'   '
	
	Do ii=0, nGridYSize-1
		WRITE(*, 200, advance='no')	'--'
	END DO
	
	WRITE(*, *)				! new line
		
100	FORMAT(I2, X)
200	FORMAT(A2, X)

	DO ii=0, nGridXSize-1
		WRITE(*, 300, advance='no')	ii, '|'
		
		DO jj=0, nGridYSize-1
			IF (gol_buf(ii, jj) == 0) THEN
				WRITE(*, 400, advance='no')	' - '
			ELSE IF (gol_buf(ii, jj) == 1) THEN
				WRITE(*, 400, advance='no')	' # '
			ELSE
				WRITE(*, 400, advance='no') ' $ '
			END IF
		END DO
		WRITE(*, *)				! new line
	END DO 

300	FORMAT(I2, A)
400	FORMAT(A)
END SUBROUTINE

FUNCTION getNeighborCount(ii, jj, tmp_buf, nGridXSize, nGridYSize)
INTEGER ii, jj, nGridXSize, nGridYSize
INTEGER, DIMENSION(0:nGridXSize, 0:nGridYSize)		:: tmp_buf
INTEGER getNeighborCount
	INTEGER a, b, c, d, e, f, x
	INTEGER count

	a = jj - 1		
	b = jj			
	c = jj + 1		
	d = ii - 1		
	e = ii			
	f = ii + 1		
	count = 0
	
	DO x=0, 7
		SELECT CASE(x)
			CASE (0)
				IF (a >=0 .AND. d >= 0)			count = count + tmp_buf(a, d)
			CASE (1)
				IF (a >= 0)				count = count + tmp_buf(a, e)
			CASE (2)
				IF (a >= 0 .AND. f < nGridYSize)	count = count + tmp_buf(a, f)
			CASE (3)
				IF (f < nGridYSize)			count = count + tmp_buf(b, f)
			CASE (4)
				IF (c < nGridXSize .AND. f < nGridYSize)count = count + tmp_buf(c, f)
			CASE (5)
				IF (c < nGridXSize)			count = count + tmp_buf(c, e)
			CASE (6)
				IF (c < nGridXSize .AND. d >= 0)	count = count + tmp_buf(c, d)
			CASE (7)
				IF (d >= 0)				count = count + tmp_buf(b, d)
		END SELECT
	END DO	

	getNeighborCount = count
END FUNCTION

SUBROUTINE calEvolution(gol_buf, nGridXSize, nGridYSize, nMyStartX, nMyStartY, nMyWidthX, nMyWidthY)
INTEGER nGridXSize, nGridYSize, nMyStartX, nMyStartY, nMyWidthX, nMyWidthY
INTEGER, DIMENSION(0:nGridXSize, 0:nGridYSize)			:: gol_buf
	INTEGER	ii, jj
	INTEGER, DIMENSION(0:nGridXSize, 0:nGridYSize)		:: tmp_buf
	INTEGER nNeighbors
	
	! Copy from gol-buf to tmp_buf
	! DO ii=0, nGridXSize-1
	!	DO jj=0, nGridYSize-1
	!		tmp_buf(ii, jj) = gol_buf(ii, jj)
	!	END DO
	! END DO

	! Copy from gol-buf to tmp_buf
	tmp_buf(0:nGridXSize, 0:nGridYSize) = gol_buf(0:nGridXSize, 0:nGridYSize)
	
	DO ii=nMyStartX, nMyStartX+nMyWidthX-1
		DO jj=nMyStartY, nMyStartY+nMyWidthY-1
			nNeighbors = getNeighborCount(jj, ii, tmp_buf, nGridXSize, nGridYSize)
			
			IF (nNeighbors == 2) THEN		! no change
				! gol_buf(ii, jj) = gol_buf(ii, jj)
			ELSE IF (nNeighbors == 3) THEN		! alive
				gol_buf(ii, jj) = 1
			ELSE					! dead
				gol_buf(ii, jj) = 0
			END IF
		END DO
	END DO
END SUBROUTINE

SUBROUTINE init_random_seed()
            INTEGER :: i, n, clock
            INTEGER, DIMENSION(:), ALLOCATABLE :: seed
          
            CALL RANDOM_SEED(size = n)
            ALLOCATE(seed(n))
          
            CALL SYSTEM_CLOCK(COUNT=clock)
          
            seed = clock + 37 * (/ (i - 1, i = 1, n) /)
            CALL RANDOM_SEED(PUT = seed)
          
            DEALLOCATE(seed)
END SUBROUTINE

!------------------------------------------------------------
END
